Summary:

Two avenues of analysis:

First we see how \(\Delta\)load affects \(\Delta\)ArcLength^2 and \(\Delta\)frequency^2, accounting for bee size and treatment order
We check to see how \(\Delta\)frequency and \(\Delta\)ArcLength are associated with \(\Delta\)MetabolicRate, while accounting for \(\Delta\)load

Measured variables:

Bee ID: Name given to an invididual (subscript in equations)
Order: This is the trial number per individual, either 1 or 2
Treatment: Either light or heavy- H is with both nectar and external load, L is just nectar load (subscript in equations)
MStarved: The empty mass of the bee after being starved for days, until it no longer buzzes when prodded in grams
M2: Total mass of the bee and load at beginning of trial in grams
MF: The mass of the bee and load at end of trial in grams
ITspan=intertegular span in m
S=area of both forewings in \(m^2\)
MetR= Metabolic rate in mL \(CO_2\) \(hr^{-1}\)
Freq= wing beat frequency in Hz
Amp= stroke amplitude of forewing in degrees
L= length of forewing from wing base to tip in m

Calculated variables:

MT: Mtrial, total mass, =(M2+MF)/2 = the average mass of the course of the trial
Load= Mtrial- Mstarved
Perload= Percent load, %load (Load/Mstarved)x100
ArcL= (.75 \(\cdotp\) wlength) (Amp (\(\pi\)/180)) pi/180 converts to radians
U= wing velocity, arcL\(\cdotp\) frq \(\cdotp\) 2 (the 2 comes from the fact that frequency consists of both an up and down stroke, vel=∆x/∆t, ∆x=arclength, ∆t=1/(freg \(\cdotp\) 2)
Freqsq=\(Freq^2\)
ArcLsq=\(ArcL^2\)
Force=\(U^2 \cdotp S\)

Install packages and read data

ipak <- function(pkg){
     new.pkg <- pkg[!(pkg %in% installed.packages()[, "Package"])]
     if(length(new.pkg)) install.packages(new.pkg, dependencies = TRUE)
     sapply(pkg, require, character.only = TRUE)
}

packages <- c('ggplot2', 'car', 'lme4', 'gsheet', "MASS",'reshape2', 
              'influence.ME', 'sjPlot', "effects", 'visreg', 'viridis')
ipak(packages)

##      ggplot2          car         lme4       gsheet         MASS 
##         TRUE         TRUE         TRUE         TRUE         TRUE 
##     reshape2 influence.ME       sjPlot      effects       visreg 
##         TRUE         TRUE         TRUE         TRUE         TRUE 
##      viridis 
##         TRUE

# read in data -- google sheet called "Final Bee Resp Data"
url = 'https://docs.google.com/spreadsheets/d/1wT-QxSJJElhhJcIXlg2hDpFKHbLyFGuqNNj2iYvA8Vo/edit?usp=sharing'
bdta <- data.frame(gsheet2tbl(url))
summary(bdta)

##     BeeID               order      Treatment            Mstarved     
##  Length:60          Min.   :1.0   Length:60          Min.   :0.0767  
##  Class :character   1st Qu.:1.0   Class :character   1st Qu.:0.1213  
##  Mode  :character   Median :1.5   Mode  :character   Median :0.1378  
##                     Mean   :1.5                      Mean   :0.1411  
##                     3rd Qu.:2.0                      3rd Qu.:0.1670  
##                     Max.   :2.0                      Max.   :0.1909  
##        M2               MF             Itspan               S            
##  Min.   :0.1062   Min.   :0.0969   Min.   :0.003860   Min.   :4.010e-05  
##  1st Qu.:0.1808   1st Qu.:0.1762   1st Qu.:0.004400   1st Qu.:5.590e-05  
##  Median :0.2191   Median :0.2107   Median :0.004620   Median :6.130e-05  
##  Mean   :0.2245   Mean   :0.2154   Mean   :0.004594   Mean   :6.141e-05  
##  3rd Qu.:0.2646   3rd Qu.:0.2576   3rd Qu.:0.004870   3rd Qu.:6.750e-05  
##  Max.   :0.3550   Max.   :0.3440   Max.   :0.005180   Max.   :7.810e-05  
##       MetR             freq            amp               L           
##  Min.   : 5.507   Min.   :163.9   Min.   : 96.88   Min.   :0.008691  
##  1st Qu.: 8.935   1st Qu.:173.0   1st Qu.:114.16   1st Qu.:0.010570  
##  Median :10.666   Median :179.1   Median :125.96   Median :0.011017  
##  Mean   :10.600   Mean   :178.9   Mean   :123.31   Mean   :0.010911  
##  3rd Qu.:12.297   3rd Qu.:185.3   3rd Qu.:133.61   3rd Qu.:0.011444  
##  Max.   :16.475   Max.   :194.7   Max.   :144.78   Max.   :0.012118

Calculate variables

bdta <- within(bdta, {
     MT = (M2 + MF)/2
     load = MT - Mstarved
     perLoad = load / Mstarved * 100
     arcL = (0.75 * L) * (amp * (pi / 180))
     U = arcL * freq * 2
     freq2 = freq^2
     arcL2 = arcL^2
     frce = U^2 * S
     
     # convert to factor variables
     order = as.factor(as.character(order))
     Treatment = as.factor(as.character(Treatment))
     BeeID = as.factor(as.character(BeeID))
})

Create a new dataframe that calculates the changes for each bee

# convert from long to wide
newDF <- data.frame()
colsTocalc = c("order", "M2", "MF", "MetR", "freq", "amp", "load", "MT", "perLoad", "frce", "arcL2", "freq2", "U", "arcL")
for(varb in colsTocalc){
     data_wide <- dcast(bdta, BeeID + Mstarved + S + Itspan + L ~ 
                             Treatment, value.var=c(varb))
     colnames(data_wide)[6:7] = paste(varb, colnames(data_wide)[6:7], sep = "_")
     if(varb == colsTocalc[1]){
       newDF <- data_wide   
     }
     else newDF <- merge(newDF, data_wide, all.y = TRUE)    
}

head(newDF)

##   BeeID Mstarved        S  Itspan          L order_H order_L   M2_H   M2_L
## 1   E01   0.1577 6.46e-05 0.00464 0.01116843       2       1 0.2634 0.1895
## 2   E03   0.1732 6.92e-05 0.00487 0.01144364       1       2 0.3106 0.2105
## 3   E05   0.1755 7.06e-05 0.00497 0.01168730       2       1 0.2776 0.2006
## 4   E09   0.1135 4.88e-05 0.00440 0.00986884       1       2 0.2282 0.1425
## 5   E10   0.1269 5.91e-05 0.00462 0.01060036       1       2 0.2166 0.1460
## 6   E11   0.1350 5.88e-05 0.00464 0.01070018       1       2 0.2223 0.1599
##     MF_H   MF_L   MetR_H    MetR_L   freq_H   freq_L    amp_H     amp_L
## 1 0.2559 0.1812 12.54241 10.393483 186.9973 179.3308 128.7145 111.49875
## 2 0.3079 0.2023 11.78387  9.214622 174.9374 167.8412 144.7784 119.98953
## 3 0.2742 0.1975 12.20102  8.981412 178.8896 166.2385 125.6778 116.22053
## 4 0.2084 0.1390  9.47114  7.581104 191.8348 190.3535 139.3246 103.30703
## 5 0.2099 0.1422 11.03371  7.159840 186.1762 169.5451 128.8257  96.88003
## 6 0.2179 0.1579 10.23580  6.806696 183.0259 168.2223 128.7916 116.18262
##    load_H  load_L    MT_H    MT_L perLoad_H perLoad_L      frce_H
## 1 0.10195 0.02765 0.25965 0.18535  64.64807  17.53329 0.003199482
## 2 0.13605 0.03320 0.30925 0.20640  78.55081  19.16859 0.003984230
## 3 0.10040 0.02355 0.27590 0.19905  57.20798  13.41880 0.003340857
## 4 0.10480 0.02725 0.21830 0.14075  92.33480  24.00881 0.002327019
## 5 0.08635 0.01720 0.21325 0.14410  68.04571  13.55398 0.002618299
## 6 0.08510 0.02390 0.22010 0.15890  63.03704  17.70370 0.002563877
##        frce_L      arcL2_H      arcL2_L  freq2_H  freq2_L      U_H
## 1 0.002208025 0.0003540922 0.0002657061 34967.99 32159.55 7.037583
## 2 0.002519159 0.0004703414 0.0003230668 30603.08 28170.67 7.587858
## 3 0.002467172 0.0003696773 0.0003161342 32001.50 27635.24 6.879020
## 4 0.001259710 0.0003239404 0.0001781022 36800.61 36234.44 6.905419
## 5 0.001228018 0.0003195386 0.0001807120 34661.58 28745.53 6.656039
## 6 0.001762569 0.0003254128 0.0002648146 33498.49 28298.73 6.603283
##        U_L     arcL_H     arcL_L
## 1 5.846363 0.01881734 0.01630050
## 2 6.033576 0.02168736 0.01797406
## 3 5.911496 0.01922699 0.01778016
## 4 5.080722 0.01799834 0.01334549
## 5 4.558361 0.01787564 0.01344292
## 6 5.475004 0.01803920 0.01627313

Calculate \(\Delta\) variables

newDF <- within(newDF, {
     deltaPercLoad = perLoad_H - perLoad_L
     avgPercLoad = (perLoad_H + perLoad_L) / 2
     deltaMetR = MetR_H - MetR_L
     deltaFrq2 = freq2_H - freq2_L
     deltaArcL = arcL_H - arcL_L
     deltaArcL2 = arcL2_H - arcL2_L
     deltaFreq2Perc = deltaFrq2 / deltaPercLoad
     deltaLoad = scale(load_H - load_L, center = TRUE, scale = FALSE)
     dLoad_nonCent = load_H - load_L
     deltaLoad2 <- deltaLoad^2
})


plot(newDF$deltaArcL2 ~ newDF$deltaFrq2)

# compare with Susie's calculations



b2 <- read.csv("~/Desktop/b2.csv")
b2 <- b2[!is.na(b2$deltaperload), ]
b2 <- b2[order(b2$BeeID, b2$Treatment), ]

plot(b2$AverageperLoad, newDF$avgPercLoad)
abline(0,1)

plot(b2$deltaload, newDF$deltaLoad)
abline(0,1)

Use PCA to combine bee sizes into a single predictor

# principle components
aa = prcomp(newDF[, c("Mstarved", "Itspan", "S", "L")], center = TRUE, scale = TRUE)
summary(aa) # 1st pc explains ~95% of the variance in the three measurements of size

## Importance of components:
##                           PC1     PC2     PC3     PC4
## Standard deviation     1.9422 0.37234 0.22746 0.19332
## Proportion of Variance 0.9431 0.03466 0.01293 0.00934
## Cumulative Proportion  0.9431 0.97772 0.99066 1.00000

biplot(aa) # shows that all three size measurement are correlated

# note, I changed the signs of the predictions so that higher PC1 values 
# correspond to bigger bees
p1 = -predict(aa)[,1] 

# add PC1 scores to dataset
newDF$size_pc1 = p1
newDF$size_pc1_2 = newDF$size_pc1^2

# show scatterplot matrix to see correlations among size predictors
car::scatterplotMatrix(newDF[, c("Mstarved", "Itspan", "S",  "L",  "size_pc1")])

1. First we see how \(\Delta\)load affects \(\Delta\)ArcLength accounting for bee size and treatment order

This is the model selection procedure that was used: 1. Fit a large model with all two-way interactions and squared terms 2. Remove non-significant predictors, starting with interactions, squared terms, and then main effects

# reformat order so that it is more interpretable
library(plyr)
newDF$order_1 <- mapvalues(newDF$order_H, from = c(2, 1), to = c("loadedSecond", "loadedFirst"))

# fit full model
m1 <- lm(deltaArcL2 ~  (deltaLoad +  size_pc1 + order_1)^2  + 
              size_pc1_2  + deltaLoad2, data = newDF)
summary(m1)

## 
## Call:
## lm(formula = deltaArcL2 ~ (deltaLoad + size_pc1 + order_1)^2 + 
##     size_pc1_2 + deltaLoad2, data = newDF)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -6.326e-05 -1.951e-05  3.244e-06  1.635e-05  4.401e-05 
## 
## Coefficients:
##                                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                    9.818e-05  1.300e-05   7.551 2.05e-07 ***
## deltaLoad                      1.859e-03  1.178e-03   1.578    0.130    
## size_pc1                      -1.804e-05  1.263e-05  -1.428    0.168    
## order_1loadedSecond            7.310e-07  1.978e-05   0.037    0.971    
## size_pc1_2                    -3.660e-06  4.905e-06  -0.746    0.464    
## deltaLoad2                     3.807e-02  4.366e-02   0.872    0.393    
## deltaLoad:size_pc1            -2.300e-06  8.387e-04  -0.003    0.998    
## deltaLoad:order_1loadedSecond -6.487e-04  2.225e-03  -0.292    0.774    
## size_pc1:order_1loadedSecond   1.204e-05  2.099e-05   0.574    0.572    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.03e-05 on 21 degrees of freedom
## Multiple R-squared:  0.575,  Adjusted R-squared:  0.413 
## F-statistic: 3.551 on 8 and 21 DF,  p-value: 0.009387

m2 <- update(m1, .~. - deltaLoad:size_pc1)
anova(m1, m2)

## Analysis of Variance Table
## 
## Model 1: deltaArcL2 ~ (deltaLoad + size_pc1 + order_1)^2 + size_pc1_2 + 
##     deltaLoad2
## Model 2: deltaArcL2 ~ deltaLoad + size_pc1 + order_1 + size_pc1_2 + deltaLoad2 + 
##     deltaLoad:order_1 + size_pc1:order_1
##   Res.Df        RSS Df   Sum of Sq  F Pr(>F)
## 1     21 1.9281e-08                         
## 2     22 1.9281e-08 -1 -6.9071e-15  0 0.9978

summary(m2)

## 
## Call:
## lm(formula = deltaArcL2 ~ deltaLoad + size_pc1 + order_1 + size_pc1_2 + 
##     deltaLoad2 + deltaLoad:order_1 + size_pc1:order_1, data = newDF)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -6.327e-05 -1.953e-05  3.263e-06  1.634e-05  4.402e-05 
## 
## Coefficients:
##                                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                    9.817e-05  1.252e-05   7.840 8.24e-08 ***
## deltaLoad                      1.860e-03  9.935e-04   1.872   0.0745 .  
## size_pc1                      -1.806e-05  7.886e-06  -2.290   0.0320 *  
## order_1loadedSecond            7.471e-07  1.846e-05   0.040   0.9681    
## size_pc1_2                    -3.672e-06  1.745e-06  -2.104   0.0470 *  
## deltaLoad2                     3.800e-02  3.336e-02   1.139   0.2669    
## deltaLoad:order_1loadedSecond -6.509e-04  2.022e-03  -0.322   0.7505    
## size_pc1:order_1loadedSecond   1.209e-05  1.187e-05   1.019   0.3194    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.96e-05 on 22 degrees of freedom
## Multiple R-squared:  0.575,  Adjusted R-squared:  0.4397 
## F-statistic: 4.251 on 7 and 22 DF,  p-value: 0.004168

m3 <- update(m2, .~. - deltaLoad:order_1)

anova(m2, m3)

## Analysis of Variance Table
## 
## Model 1: deltaArcL2 ~ deltaLoad + size_pc1 + order_1 + size_pc1_2 + deltaLoad2 + 
##     deltaLoad:order_1 + size_pc1:order_1
## Model 2: deltaArcL2 ~ deltaLoad + size_pc1 + order_1 + size_pc1_2 + deltaLoad2 + 
##     size_pc1:order_1
##   Res.Df        RSS Df   Sum of Sq      F Pr(>F)
## 1     22 1.9281e-08                             
## 2     23 1.9371e-08 -1 -9.0841e-11 0.1037 0.7505

summary(m3)

## 
## Call:
## lm(formula = deltaArcL2 ~ deltaLoad + size_pc1 + order_1 + size_pc1_2 + 
##     deltaLoad2 + size_pc1:order_1, data = newDF)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -6.414e-05 -1.856e-05  3.421e-06  1.715e-05  4.272e-05 
## 
## Coefficients:
##                                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                   9.768e-05  1.218e-05   8.017 4.13e-08 ***
## deltaLoad                     1.625e-03  6.585e-04   2.467   0.0215 *  
## size_pc1                     -1.728e-05  7.355e-06  -2.349   0.0278 *  
## order_1loadedSecond           3.683e-06  1.573e-05   0.234   0.8170    
## size_pc1_2                   -3.816e-06  1.654e-06  -2.308   0.0304 *  
## deltaLoad2                    4.504e-02  2.467e-02   1.825   0.0809 .  
## size_pc1:order_1loadedSecond  1.003e-05  9.807e-06   1.023   0.3169    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.902e-05 on 23 degrees of freedom
## Multiple R-squared:  0.573,  Adjusted R-squared:  0.4615 
## F-statistic: 5.143 on 6 and 23 DF,  p-value: 0.001756

m4 <- update(m3, .~. - size_pc1:order_1)
anova(m3, m4)

## Analysis of Variance Table
## 
## Model 1: deltaArcL2 ~ deltaLoad + size_pc1 + order_1 + size_pc1_2 + deltaLoad2 + 
##     size_pc1:order_1
## Model 2: deltaArcL2 ~ deltaLoad + size_pc1 + order_1 + size_pc1_2 + deltaLoad2
##   Res.Df        RSS Df   Sum of Sq      F Pr(>F)
## 1     23 1.9371e-08                             
## 2     24 2.0253e-08 -1 -8.8175e-10 1.0469 0.3169

summary(m4)

## 
## Call:
## lm(formula = deltaArcL2 ~ deltaLoad + size_pc1 + order_1 + size_pc1_2 + 
##     deltaLoad2, data = newDF)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -6.427e-05 -1.874e-05 -5.000e-09  1.783e-05  4.254e-05 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)          1.024e-04  1.128e-05   9.079 3.13e-09 ***
## deltaLoad            1.639e-03  6.590e-04   2.486   0.0203 *  
## size_pc1            -1.180e-05  5.042e-06  -2.340   0.0279 *  
## order_1loadedSecond -5.411e-07  1.519e-05  -0.036   0.9719    
## size_pc1_2          -2.627e-06  1.178e-06  -2.231   0.0353 *  
## deltaLoad2           2.613e-02  1.636e-02   1.597   0.1233    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.905e-05 on 24 degrees of freedom
## Multiple R-squared:  0.5535, Adjusted R-squared:  0.4605 
## F-statistic: 5.951 on 5 and 24 DF,  p-value: 0.001027

m5 <- update(m4, .~. - deltaLoad2)
anova(m4,m5)

## Analysis of Variance Table
## 
## Model 1: deltaArcL2 ~ deltaLoad + size_pc1 + order_1 + size_pc1_2 + deltaLoad2
## Model 2: deltaArcL2 ~ deltaLoad + size_pc1 + order_1 + size_pc1_2
##   Res.Df        RSS Df   Sum of Sq      F Pr(>F)
## 1     24 2.0253e-08                             
## 2     25 2.2406e-08 -1 -2.1529e-09 2.5511 0.1233

summary(m5)

## 
## Call:
## lm(formula = deltaArcL2 ~ deltaLoad + size_pc1 + order_1 + size_pc1_2, 
##     data = newDF)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -6.782e-05 -1.748e-05  9.200e-07  1.678e-05  4.278e-05 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)          1.074e-04  1.117e-05   9.623 6.93e-10 ***
## deltaLoad            2.095e-03  6.122e-04   3.422  0.00215 ** 
## size_pc1            -1.388e-05  5.018e-06  -2.767  0.01050 *  
## order_1loadedSecond  2.783e-07  1.565e-05   0.018  0.98595    
## size_pc1_2          -2.037e-06  1.152e-06  -1.768  0.08930 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.994e-05 on 25 degrees of freedom
## Multiple R-squared:  0.5061, Adjusted R-squared:  0.427 
## F-statistic: 6.403 on 4 and 25 DF,  p-value: 0.001086

m6 <- update(m5, .~.  - size_pc1_2)
anova(m6, m5)

## Analysis of Variance Table
## 
## Model 1: deltaArcL2 ~ deltaLoad + size_pc1 + order_1
## Model 2: deltaArcL2 ~ deltaLoad + size_pc1 + order_1 + size_pc1_2
##   Res.Df        RSS Df  Sum of Sq      F Pr(>F)  
## 1     26 2.5207e-08                              
## 2     25 2.2406e-08  1 2.8009e-09 3.1251 0.0893 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

summary(m6)

## 
## Call:
## lm(formula = deltaArcL2 ~ deltaLoad + size_pc1 + order_1, data = newDF)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -6.154e-05 -1.781e-05 -4.100e-07  1.849e-05  4.983e-05 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)          9.844e-05  1.034e-05   9.525 5.79e-10 ***
## deltaLoad            2.148e-03  6.359e-04   3.378  0.00231 ** 
## size_pc1            -1.321e-05  5.205e-06  -2.539  0.01745 *  
## order_1loadedSecond  3.224e-06  1.618e-05   0.199  0.84364    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.114e-05 on 26 degrees of freedom
## Multiple R-squared:  0.4443, Adjusted R-squared:  0.3802 
## F-statistic: 6.929 on 3 and 26 DF,  p-value: 0.001402

m7 <- update(m6, .~. - order_1)
anova(m7, m6)

## Analysis of Variance Table
## 
## Model 1: deltaArcL2 ~ deltaLoad + size_pc1
## Model 2: deltaArcL2 ~ deltaLoad + size_pc1 + order_1
##   Res.Df        RSS Df  Sum of Sq      F Pr(>F)
## 1     27 2.5245e-08                            
## 2     26 2.5207e-08  1 3.8477e-11 0.0397 0.8436

summary(m7)

## 
## Call:
## lm(formula = deltaArcL2 ~ deltaLoad + size_pc1, data = newDF)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -6.041e-05 -1.877e-05 -2.970e-07  1.831e-05  4.778e-05 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.002e-04  5.583e-06  17.940  < 2e-16 ***
## deltaLoad    2.059e-03  4.439e-04   4.638 8.06e-05 ***
## size_pc1    -1.256e-05  3.967e-06  -3.166  0.00381 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.058e-05 on 27 degrees of freedom
## Multiple R-squared:  0.4435, Adjusted R-squared:  0.4022 
## F-statistic: 10.76 on 2 and 27 DF,  p-value: 0.0003666

anova(m7, update(m7, .~. - deltaLoad)) # p-value for delta load

## Analysis of Variance Table
## 
## Model 1: deltaArcL2 ~ deltaLoad + size_pc1
## Model 2: deltaArcL2 ~ size_pc1
##   Res.Df        RSS Df   Sum of Sq      F    Pr(>F)    
## 1     27 2.5245e-08                                    
## 2     28 4.5359e-08 -1 -2.0114e-08 21.512 8.055e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

anova(m7, update(m7, .~. - size_pc1)) # p-value for size

## Analysis of Variance Table
## 
## Model 1: deltaArcL2 ~ deltaLoad + size_pc1
## Model 2: deltaArcL2 ~ deltaLoad
##   Res.Df        RSS Df   Sum of Sq      F   Pr(>F)   
## 1     27 2.5245e-08                                  
## 2     28 3.4619e-08 -1 -9.3737e-09 10.025 0.003808 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

summary(m7)

## 
## Call:
## lm(formula = deltaArcL2 ~ deltaLoad + size_pc1, data = newDF)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -6.041e-05 -1.877e-05 -2.970e-07  1.831e-05  4.778e-05 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.002e-04  5.583e-06  17.940  < 2e-16 ***
## deltaLoad    2.059e-03  4.439e-04   4.638 8.06e-05 ***
## size_pc1    -1.256e-05  3.967e-06  -3.166  0.00381 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.058e-05 on 27 degrees of freedom
## Multiple R-squared:  0.4435, Adjusted R-squared:  0.4022 
## F-statistic: 10.76 on 2 and 27 DF,  p-value: 0.0003666

# refit model with non-centered version of deltaLoad
# this model will have a different intercept, but same p-values for slopes
summary(update(m7, .~. - deltaLoad + dLoad_nonCent)) # final model for paper

## 
## Call:
## lm(formula = deltaArcL2 ~ size_pc1 + dLoad_nonCent, data = newDF)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -6.041e-05 -1.877e-05 -2.970e-07  1.831e-05  4.778e-05 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   -5.822e-05  3.460e-05  -1.683  0.10397    
## size_pc1      -1.256e-05  3.967e-06  -3.166  0.00381 ** 
## dLoad_nonCent  2.059e-03  4.439e-04   4.638 8.06e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.058e-05 on 27 degrees of freedom
## Multiple R-squared:  0.4435, Adjusted R-squared:  0.4022 
## F-statistic: 10.76 on 2 and 27 DF,  p-value: 0.0003666

model diagnostics

par(mfrow = c(2,2))
plot(m7, which = 1:4) # no glaring violations, though there are a few fairly influential observations

par(mfrow = c(1,1))

car::vif(m7)

## deltaLoad  size_pc1 
##  1.841066  1.841066

Model visualization

summary(m7)

## 
## Call:
## lm(formula = deltaArcL2 ~ deltaLoad + size_pc1, data = newDF)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -6.041e-05 -1.877e-05 -2.970e-07  1.831e-05  4.778e-05 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.002e-04  5.583e-06  17.940  < 2e-16 ***
## deltaLoad    2.059e-03  4.439e-04   4.638 8.06e-05 ***
## size_pc1    -1.256e-05  3.967e-06  -3.166  0.00381 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.058e-05 on 27 degrees of freedom
## Multiple R-squared:  0.4435, Adjusted R-squared:  0.4022 
## F-statistic: 10.76 on 2 and 27 DF,  p-value: 0.0003666

# calculate partial residuals for deltaLoad
# these are the  residuals, minus the effect of detlaLoad
y <- residuals(m7, type = 'partial')[, "deltaLoad"]

# plot partial residuals with base R plotting
plot(x = newDF$deltaLoad, y = y)

# double check to make sure the slope for partial residual plots are the 
# same as in the original regression
summary(lm(y ~ newDF$dLoad_nonCent))

## 
## Call:
## lm(formula = y ~ newDF$dLoad_nonCent)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -6.041e-05 -1.877e-05 -2.970e-07  1.831e-05  4.778e-05 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)         -1.584e-04  2.531e-05  -6.257 9.20e-07 ***
## newDF$dLoad_nonCent  2.059e-03  3.213e-04   6.409 6.14e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.003e-05 on 28 degrees of freedom
## Multiple R-squared:  0.5946, Adjusted R-squared:  0.5801 
## F-statistic: 41.07 on 1 and 28 DF,  p-value: 6.137e-07

# this is what the raw data look like
plot(x = newDF$deltaLoad, y = newDF$deltaFrq2)

# this package plots the partial residuals
crPlot(m7, variable = "deltaLoad")

# plot with ggplot2
theme_set(theme_classic()) 

# plot raw data w/ ggplot
ggplot(newDF, aes(x= size_pc1, y = deltaArcL, color = deltaLoad)) + 
     geom_point()

ggplot(newDF, aes(x= deltaLoad, y = deltaArcL, color = size_pc1)) + 
     geom_point()

# y axis isn't easily interpretable
ggplot(newDF, aes(x= dLoad_nonCent, y = y)) + 
     geom_point() + 
     labs(x = "delta load", y = "partial residuals for delta load \n i.e. delta load effect on arcL^2 \n while subtracing affect of bee size") + 
     stat_smooth(method = 'lm', se = FALSE)

Takeaways from model 1:

A higher deltaload (i.e. a relatively larger “high” load) causes bees to have a larger change in arclength squared. Another way to say this is that increasing the change in load by 1 gram causes an increase in the change in arclength^2 by 0.002059 degrees.
The larger the bee, the lower the change in arcLength^2 (association, rather than causation), while holding the change in load constant.
We do not have enough evidence to notice any non-linear relationships in this model.

See how \(\Delta\)load affects \(\Delta\)frequency, accounting for bee size and treatment order

# fit full model
m1 <- lm(deltaFrq2 ~  (deltaLoad +  size_pc1 + order_1)^2  + 
              size_pc1_2  + deltaLoad2, data = newDF)
summary(m1)

## 
## Call:
## lm(formula = deltaFrq2 ~ (deltaLoad + size_pc1 + order_1)^2 + 
##     size_pc1_2 + deltaLoad2, data = newDF)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3981.4 -1255.9  -170.3  1165.6  3552.9 
## 
## Coefficients:
##                                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                      4195.97     909.39   4.614  0.00015 ***
## deltaLoad                     -106368.80   82390.69  -1.291  0.21073    
## size_pc1                          666.06     883.26   0.754  0.45917    
## order_1loadedSecond             -4376.19    1383.48  -3.163  0.00469 ** 
## size_pc1_2                       -120.74     343.06  -0.352  0.72839    
## deltaLoad2                    -469149.49 3053199.27  -0.154  0.87935    
## deltaLoad:size_pc1               6902.45   58657.40   0.118  0.90744    
## deltaLoad:order_1loadedSecond   62840.01  155625.35   0.404  0.69045    
## size_pc1:order_1loadedSecond       54.39    1467.95   0.037  0.97079    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2119 on 21 degrees of freedom
## Multiple R-squared:  0.5048, Adjusted R-squared:  0.3162 
## F-statistic: 2.676 on 8 and 21 DF,  p-value: 0.03372

car::vif(m1) # some serious multicollinearity

##          deltaLoad           size_pc1            order_1 
##          13.203325          19.003708           3.182227 
##         size_pc1_2         deltaLoad2 deltaLoad:size_pc1 
##          18.740891          10.117446          25.115087 
##  deltaLoad:order_1   size_pc1:order_1 
##          11.384208          23.301468

m2 <- update(m1, .~. - size_pc1:order_1)
anova(m1, m2)

## Analysis of Variance Table
## 
## Model 1: deltaFrq2 ~ (deltaLoad + size_pc1 + order_1)^2 + size_pc1_2 + 
##     deltaLoad2
## Model 2: deltaFrq2 ~ deltaLoad + size_pc1 + order_1 + size_pc1_2 + deltaLoad2 + 
##     deltaLoad:size_pc1 + deltaLoad:order_1
##   Res.Df      RSS Df Sum of Sq      F Pr(>F)
## 1     21 94310333                           
## 2     22 94316499 -1   -6165.2 0.0014 0.9708

summary(m2)

## 
## Call:
## lm(formula = deltaFrq2 ~ deltaLoad + size_pc1 + order_1 + size_pc1_2 + 
##     deltaLoad2 + deltaLoad:size_pc1 + deltaLoad:order_1, data = newDF)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3973.1 -1259.6  -176.2  1149.9  3556.7 
## 
## Coefficients:
##                                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                      4205.4      852.7   4.932 6.21e-05 ***
## deltaLoad                     -108212.1    64166.1  -1.686  0.10584    
## size_pc1                          695.7      367.7   1.892  0.07176 .  
## order_1loadedSecond             -4386.5     1324.0  -3.313  0.00316 ** 
## size_pc1_2                       -109.9      176.0  -0.625  0.53864    
## deltaLoad2                    -417980.4  2660489.9  -0.157  0.87659    
## deltaLoad:size_pc1               5130.1    33170.6   0.155  0.87850    
## deltaLoad:order_1loadedSecond   66238.1   122843.3   0.539  0.59516    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2071 on 22 degrees of freedom
## Multiple R-squared:  0.5048, Adjusted R-squared:  0.3472 
## F-statistic: 3.203 on 7 and 22 DF,  p-value: 0.01701

m3 <- update(m2, .~. - deltaLoad:size_pc1)

anova(m2, m3)

## Analysis of Variance Table
## 
## Model 1: deltaFrq2 ~ deltaLoad + size_pc1 + order_1 + size_pc1_2 + deltaLoad2 + 
##     deltaLoad:size_pc1 + deltaLoad:order_1
## Model 2: deltaFrq2 ~ deltaLoad + size_pc1 + order_1 + size_pc1_2 + deltaLoad2 + 
##     deltaLoad:order_1
##   Res.Df      RSS Df Sum of Sq      F Pr(>F)
## 1     22 94316499                           
## 2     23 94419044 -1   -102545 0.0239 0.8785

summary(m3)

## 
## Call:
## lm(formula = deltaFrq2 ~ deltaLoad + size_pc1 + order_1 + size_pc1_2 + 
##     deltaLoad2 + deltaLoad:order_1, data = newDF)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3910.3 -1244.0  -168.5  1173.4  3503.8 
## 
## Coefficients:
##                                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                      4190.99     829.34   5.053 4.09e-05 ***
## deltaLoad                     -108478.92   62767.13  -1.728  0.09734 .  
## size_pc1                          702.89     356.93   1.969  0.06108 .  
## order_1loadedSecond             -4434.40    1259.70  -3.520  0.00184 ** 
## size_pc1_2                        -88.67     107.58  -0.824  0.41829    
## deltaLoad2                    -215246.26 2265485.45  -0.095  0.92513    
## deltaLoad:order_1loadedSecond   61643.10  116639.62   0.528  0.60222    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2026 on 23 degrees of freedom
## Multiple R-squared:  0.5042, Adjusted R-squared:  0.3749 
## F-statistic: 3.899 on 6 and 23 DF,  p-value: 0.00785

m4 <- update(m3, .~. - deltaLoad:order_1)
anova(m3, m4)

## Analysis of Variance Table
## 
## Model 1: deltaFrq2 ~ deltaLoad + size_pc1 + order_1 + size_pc1_2 + deltaLoad2 + 
##     deltaLoad:order_1
## Model 2: deltaFrq2 ~ deltaLoad + size_pc1 + order_1 + size_pc1_2 + deltaLoad2
##   Res.Df      RSS Df Sum of Sq      F Pr(>F)
## 1     23 94419044                           
## 2     24 95565635 -1  -1146591 0.2793 0.6022

summary(m4)

## 
## Call:
## lm(formula = deltaFrq2 ~ deltaLoad + size_pc1 + order_1 + size_pc1_2 + 
##     deltaLoad2, data = newDF)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3667.1 -1308.7   -52.9  1317.9  3482.6 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)          4.330e+03  7.749e+02   5.587 9.48e-06 ***
## deltaLoad           -8.589e+04  4.527e+04  -1.897 0.069889 .  
## size_pc1             7.352e+02  3.463e+02   2.123 0.044282 *  
## order_1loadedSecond -4.794e+03  1.044e+03  -4.594 0.000117 ***
## size_pc1_2          -5.195e+01  8.090e+01  -0.642 0.526837    
## deltaLoad2          -1.250e+06  1.124e+06  -1.112 0.277131    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1995 on 24 degrees of freedom
## Multiple R-squared:  0.4982, Adjusted R-squared:  0.3937 
## F-statistic: 4.766 on 5 and 24 DF,  p-value: 0.00364

m5 <- update(m4, .~. - size_pc1_2)
anova(m4,m5)

## Analysis of Variance Table
## 
## Model 1: deltaFrq2 ~ deltaLoad + size_pc1 + order_1 + size_pc1_2 + deltaLoad2
## Model 2: deltaFrq2 ~ deltaLoad + size_pc1 + order_1 + deltaLoad2
##   Res.Df      RSS Df Sum of Sq      F Pr(>F)
## 1     24 95565635                           
## 2     25 97207799 -1  -1642165 0.4124 0.5268

summary(m5)

## 
## Call:
## lm(formula = deltaFrq2 ~ deltaLoad + size_pc1 + order_1 + deltaLoad2, 
##     data = newDF)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3591.2 -1388.0  -154.6  1391.9  3336.4 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)             4166.0      723.2   5.760 5.29e-06 ***
## deltaLoad             -80714.2    44020.7  -1.834 0.078654 .  
## size_pc1                 732.5      342.2   2.141 0.042251 *  
## order_1loadedSecond    -4719.5     1024.9  -4.605 0.000104 ***
## deltaLoad2          -1475818.2  1054446.2  -1.400 0.173916    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1972 on 25 degrees of freedom
## Multiple R-squared:  0.4896, Adjusted R-squared:  0.4079 
## F-statistic: 5.995 on 4 and 25 DF,  p-value: 0.00159

m6 <- update(m5, .~.  - deltaLoad2)
anova(m6, m5)

## Analysis of Variance Table
## 
## Model 1: deltaFrq2 ~ deltaLoad + size_pc1 + order_1
## Model 2: deltaFrq2 ~ deltaLoad + size_pc1 + order_1 + deltaLoad2
##   Res.Df       RSS Df Sum of Sq      F Pr(>F)
## 1     26 104824695                           
## 2     25  97207799  1   7616896 1.9589 0.1739

summary(m6)

## 
## Call:
## lm(formula = deltaFrq2 ~ deltaLoad + size_pc1 + order_1, data = newDF)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3191.2 -1541.4  -264.1  1651.1  3486.0 
## 
## Coefficients:
##                      Estimate Std. Error t value Pr(>|t|)    
## (Intercept)            3735.3      666.5   5.605 6.87e-06 ***
## deltaLoad           -105594.7    41007.6  -2.575 0.016065 *  
## size_pc1                861.4      335.6   2.566 0.016382 *  
## order_1loadedSecond   -4717.7     1043.6  -4.521 0.000119 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2008 on 26 degrees of freedom
## Multiple R-squared:  0.4496, Adjusted R-squared:  0.3861 
## F-statistic:  7.08 on 3 and 26 DF,  p-value: 0.001244

m7 <- update(m6, .~. - size_pc1)
anova(m7, m6) # p-values for size

## Analysis of Variance Table
## 
## Model 1: deltaFrq2 ~ deltaLoad + order_1
## Model 2: deltaFrq2 ~ deltaLoad + size_pc1 + order_1
##   Res.Df       RSS Df Sum of Sq      F  Pr(>F)  
## 1     27 131380648                              
## 2     26 104824695  1  26555953 6.5868 0.01638 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

m8 <- update(m6, .~. - deltaLoad)
anova(m6, m8) # p-value for deltaLoad

## Analysis of Variance Table
## 
## Model 1: deltaFrq2 ~ deltaLoad + size_pc1 + order_1
## Model 2: deltaFrq2 ~ size_pc1 + order_1
##   Res.Df       RSS Df Sum of Sq      F  Pr(>F)  
## 1     26 104824695                              
## 2     27 131557512 -1 -26732816 6.6306 0.01607 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

m9 <- update(m6, .~. - order_1)
anova(m6, m9) # p-value for order

## Analysis of Variance Table
## 
## Model 1: deltaFrq2 ~ deltaLoad + size_pc1 + order_1
## Model 2: deltaFrq2 ~ deltaLoad + size_pc1
##   Res.Df       RSS Df Sum of Sq      F    Pr(>F)    
## 1     26 104824695                                  
## 2     27 187214161 -1 -82389465 20.435 0.0001192 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# refit final model with non-centered deltaLoad
m10 <- update(m6, .~. - deltaLoad + dLoad_nonCent)
summary(m10) # final model for paper

## 
## Call:
## lm(formula = deltaFrq2 ~ size_pc1 + order_1 + dLoad_nonCent, 
##     data = newDF)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3191.2 -1541.4  -264.1  1651.1  3486.0 
## 
## Coefficients:
##                      Estimate Std. Error t value Pr(>|t|)    
## (Intercept)           11857.7     3586.6   3.306 0.002766 ** 
## size_pc1                861.4      335.6   2.566 0.016382 *  
## order_1loadedSecond   -4717.7     1043.6  -4.521 0.000119 ***
## dLoad_nonCent       -105594.7    41007.6  -2.575 0.016065 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2008 on 26 degrees of freedom
## Multiple R-squared:  0.4496, Adjusted R-squared:  0.3861 
## F-statistic:  7.08 on 3 and 26 DF,  p-value: 0.001244

model diagnostics

par(mfrow = c(2,2))
plot(m10, which = 1:4) # no glaring violations

par(mfrow = c(1,1))

car::vif(m10) # vif is a little high

##      size_pc1       order_1 dLoad_nonCent 
##      3.056456      2.017003      3.643394

model visualization

# calculate partial residuals for deltaLoad
# these are the  residuals, minus the effect of detlaLoad
y <- residuals(m10, type = 'partial')[, "dLoad_nonCent"]

# plot partial residuals with base R plotting
plot(x = newDF$dLoad_nonCent, y = y)

# double check to make sure the slope for partial residual plots are the 
# same as in the original regression
summary(lm(y ~ newDF$dLoad_nonCent))

## 
## Call:
## lm(formula = y ~ newDF$dLoad_nonCent)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3191.2 -1541.4  -264.1  1651.1  3486.0 
## 
## Coefficients:
##                     Estimate Std. Error t value Pr(>|t|)    
## (Intercept)             8122       1631   4.980 2.93e-05 ***
## newDF$dLoad_nonCent  -105595      20702  -5.101 2.11e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1935 on 28 degrees of freedom
## Multiple R-squared:  0.4816, Adjusted R-squared:  0.4631 
## F-statistic: 26.02 on 1 and 28 DF,  p-value: 2.106e-05

# this is what the raw data look like
plot(x = newDF$dLoad_nonCent, y = newDF$deltaFrq2)

# this package plots the partial residuals
crPlot(m10, variable = "dLoad_nonCent")

# plot with ggplot2
# plot raw data w/ ggplot
ggplot(newDF, aes(x= size_pc1, y = deltaFrq2, color = deltaLoad, shape = order_1)) + 
     geom_point()

ggplot(newDF, aes(x= deltaLoad, y = deltaFrq2, color = size_pc1, shape = order_1)) + 
     geom_point()

# y axis isn't easily interpretable
ggplot(newDF, aes(x= dLoad_nonCent, y = y)) + 
     geom_point() + 
     labs(x = "delta load", y = "partial residuals for delta load \n i.e. delta load effect on freq^2 \n while subtracing affect of bee size and order") + 
     stat_smooth(method = 'lm', se = FALSE)

# partial residuals for order

y <- residuals(m10, type = 'partial')[, "order_1"]

ggplot(newDF, aes(x= order_1, y = y)) + 
     geom_boxplot() + 
     labs(x = "order", y = "partial residuals for order \n i.e. order effect on freq^2 \n while subtracing affect of bee size and load") + 
     stat_smooth(method = 'lm', se = FALSE)

# holding bee size and delta load constant, if a bee was loaded second, (confusing, huh?) then it would have a much lower freq^2 than if it was loaded first

summary(m10)

## 
## Call:
## lm(formula = deltaFrq2 ~ size_pc1 + order_1 + dLoad_nonCent, 
##     data = newDF)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3191.2 -1541.4  -264.1  1651.1  3486.0 
## 
## Coefficients:
##                      Estimate Std. Error t value Pr(>|t|)    
## (Intercept)           11857.7     3586.6   3.306 0.002766 ** 
## size_pc1                861.4      335.6   2.566 0.016382 *  
## order_1loadedSecond   -4717.7     1043.6  -4.521 0.000119 ***
## dLoad_nonCent       -105594.7    41007.6  -2.575 0.016065 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2008 on 26 degrees of freedom
## Multiple R-squared:  0.4496, Adjusted R-squared:  0.3861 
## F-statistic:  7.08 on 3 and 26 DF,  p-value: 0.001244

Takeaways from model 2:

Order, size, and deltaLoad are associated with a change in frequency^2. We find no evidence of non-linear relationships or interactions.
Holding size and order constant, we find that a larger deltaLoad causes a decrease in deltaFrequency^2
Holding other variables constant, an increase in bee size is associated with a larger deltaFrequency^2.
Holding other variables constant, if the bee was loaded second, then they had a lower deltaFrequency^2 than if they were loaded first.

# make some predictions
predDF <- data.frame(size_pc1 = 0, order_1 = c("loadedSecond", "loadedFirst"), dLoad_nonCent = mean(newDF$dLoad_nonCent))

predDF$pred_deltaF2 <- predict(m10, predDF)

predDF

##   size_pc1      order_1 dLoad_nonCent pred_deltaF2
## 1        0 loadedSecond       0.07692     -982.314
## 2        0  loadedFirst       0.07692     3735.337

2. We check to see how \(\Delta\)frequency and \(\Delta\)ArcLength are associated with \(\Delta\)MetabolicRate, while accounting for \(\Delta\)load

mm1 <- lm(deltaMetR ~ deltaFrq2 + deltaArcL2 + deltaLoad + size_pc1 + deltaLoad2, data = newDF)
car::vif(mm1)

##  deltaFrq2 deltaArcL2  deltaLoad   size_pc1 deltaLoad2 
##   1.060059   1.852143   3.844322   2.679591   1.495113

summary(mm1)

## 
## Call:
## lm(formula = deltaMetR ~ deltaFrq2 + deltaArcL2 + deltaLoad + 
##     size_pc1 + deltaLoad2, data = newDF)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.11198 -0.64060 -0.01055  0.42457  2.97948 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.941e+00  6.182e-01   3.140  0.00444 ** 
## deltaFrq2    3.543e-04  6.888e-05   5.144 2.89e-05 ***
## deltaArcL2  -2.573e+03  5.900e+03  -0.436  0.66659    
## deltaLoad    2.368e+01  1.937e+01   1.223  0.23337    
## size_pc1    -5.482e-02  1.445e-01  -0.379  0.70776    
## deltaLoad2   4.440e+02  5.114e+02   0.868  0.39381    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.9233 on 24 degrees of freedom
## Multiple R-squared:  0.5937, Adjusted R-squared:  0.5091 
## F-statistic: 7.014 on 5 and 24 DF,  p-value: 0.0003629

mm2 <- update(mm1, .~. - size_pc1)
anova(mm1, mm2)

## Analysis of Variance Table
## 
## Model 1: deltaMetR ~ deltaFrq2 + deltaArcL2 + deltaLoad + size_pc1 + deltaLoad2
## Model 2: deltaMetR ~ deltaFrq2 + deltaArcL2 + deltaLoad + deltaLoad2
##   Res.Df    RSS Df Sum of Sq      F Pr(>F)
## 1     24 20.460                           
## 2     25 20.582 -1  -0.12268 0.1439 0.7078

summary(mm2)

## 
## Call:
## lm(formula = deltaMetR ~ deltaFrq2 + deltaArcL2 + deltaLoad + 
##     deltaLoad2, data = newDF)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.06884 -0.60009 -0.03648  0.36122  3.09536 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.821e+00  5.223e-01   3.487  0.00182 ** 
## deltaFrq2    3.571e-04  6.730e-05   5.306 1.69e-05 ***
## deltaArcL2  -1.543e+03  5.147e+03  -0.300  0.76682    
## deltaLoad    1.788e+01  1.168e+01   1.530  0.13849    
## deltaLoad2   4.891e+02  4.888e+02   1.001  0.32665    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.9074 on 25 degrees of freedom
## Multiple R-squared:  0.5913, Adjusted R-squared:  0.5259 
## F-statistic: 9.041 on 4 and 25 DF,  p-value: 0.0001166

mm3 <- update(mm2, .~. - deltaArcL2)
anova(mm3, mm2)

## Analysis of Variance Table
## 
## Model 1: deltaMetR ~ deltaFrq2 + deltaLoad + deltaLoad2
## Model 2: deltaMetR ~ deltaFrq2 + deltaArcL2 + deltaLoad + deltaLoad2
##   Res.Df    RSS Df Sum of Sq      F Pr(>F)
## 1     26 20.656                           
## 2     25 20.582  1   0.07399 0.0899 0.7668

summary(mm3)

## 
## Call:
## lm(formula = deltaMetR ~ deltaFrq2 + deltaLoad + deltaLoad2, 
##     data = newDF)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.05850 -0.58135 -0.00991  0.37733  3.10956 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 1.682e+00  2.327e-01   7.228 1.12e-07 ***
## deltaFrq2   3.560e-04  6.602e-05   5.393 1.19e-05 ***
## deltaLoad   1.667e+01  1.077e+01   1.548    0.134    
## deltaLoad2  4.422e+02  4.550e+02   0.972    0.340    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.8913 on 26 degrees of freedom
## Multiple R-squared:  0.5898, Adjusted R-squared:  0.5425 
## F-statistic: 12.46 on 3 and 26 DF,  p-value: 3.065e-05

mm4 <- update(mm3, .~. - deltaLoad2)
anova(mm3, mm4)

## Analysis of Variance Table
## 
## Model 1: deltaMetR ~ deltaFrq2 + deltaLoad + deltaLoad2
## Model 2: deltaMetR ~ deltaFrq2 + deltaLoad
##   Res.Df    RSS Df Sum of Sq      F Pr(>F)
## 1     26 20.656                           
## 2     27 21.407 -1  -0.75053 0.9447   0.34

summary(mm4)

## 
## Call:
## lm(formula = deltaMetR ~ deltaFrq2 + deltaLoad, data = newDF)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.17958 -0.59123 -0.06454  0.34379  3.00903 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 1.824e+00  1.808e-01  10.086 1.18e-10 ***
## deltaFrq2   3.451e-04  6.498e-05   5.310 1.32e-05 ***
## deltaLoad   2.140e+01  9.595e+00   2.231   0.0342 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.8904 on 27 degrees of freedom
## Multiple R-squared:  0.5749, Adjusted R-squared:  0.5434 
## F-statistic: 18.26 on 2 and 27 DF,  p-value: 9.654e-06

mm5 <- update(mm4, .~. - deltaLoad) 
anova(mm5, mm4) # p-value for load

## Analysis of Variance Table
## 
## Model 1: deltaMetR ~ deltaFrq2
## Model 2: deltaMetR ~ deltaFrq2 + deltaLoad
##   Res.Df    RSS Df Sum of Sq      F  Pr(>F)  
## 1     28 25.352                              
## 2     27 21.407  1    3.9448 4.9754 0.03421 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

mm6 <- update(mm4, .~. - deltaFrq2)
anova(mm4, mm6) # p-value for deltafrq2

## Analysis of Variance Table
## 
## Model 1: deltaMetR ~ deltaFrq2 + deltaLoad
## Model 2: deltaMetR ~ deltaLoad
##   Res.Df    RSS Df Sum of Sq      F    Pr(>F)    
## 1     27 21.407                                  
## 2     28 43.766 -1   -22.359 28.201 1.324e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# refit model with non-centered load
mm7 <- update(mm4, .~. - deltaLoad + dLoad_nonCent)

summary(mm7) # final model for paper

## 
## Call:
## lm(formula = deltaMetR ~ deltaFrq2 + dLoad_nonCent, data = newDF)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.17958 -0.59123 -0.06454  0.34379  3.00903 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   1.777e-01  7.507e-01   0.237   0.8146    
## deltaFrq2     3.451e-04  6.498e-05   5.310 1.32e-05 ***
## dLoad_nonCent 2.140e+01  9.595e+00   2.231   0.0342 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.8904 on 27 degrees of freedom
## Multiple R-squared:  0.5749, Adjusted R-squared:  0.5434 
## F-statistic: 18.26 on 2 and 27 DF,  p-value: 9.654e-06

model diagnostics

car::vif(mm7)

##     deltaFrq2 dLoad_nonCent 
##      1.014369      1.014369

par(mfrow = c(2,2))
plot(mm7, which = 1:4) # no glaring violations

par(mfrow = c(1,1))

model visualization

# calculate partial residuals for deltaLoad
# these are the  residuals, minus the effect of detlaLoad
y <- residuals(mm7, type = 'partial')[, "dLoad_nonCent"]

# plot partial residuals with base R plotting
plot(x = newDF$dLoad_nonCent, y = y)

# double check to make sure the slope for partial residual plots are the 
# same as in the original regression
summary(lm(y ~ newDF$dLoad_nonCent))

## 
## Call:
## lm(formula = y ~ newDF$dLoad_nonCent)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.17958 -0.59123 -0.06454  0.34379  3.00903 
## 
## Coefficients:
##                     Estimate Std. Error t value Pr(>|t|)  
## (Intercept)          -1.6463     0.7371  -2.233   0.0337 *
## newDF$dLoad_nonCent  21.4030     9.3554   2.288   0.0299 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.8744 on 28 degrees of freedom
## Multiple R-squared:  0.1575, Adjusted R-squared:  0.1274 
## F-statistic: 5.234 on 1 and 28 DF,  p-value: 0.02991

summary(mm7)

## 
## Call:
## lm(formula = deltaMetR ~ deltaFrq2 + dLoad_nonCent, data = newDF)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.17958 -0.59123 -0.06454  0.34379  3.00903 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   1.777e-01  7.507e-01   0.237   0.8146    
## deltaFrq2     3.451e-04  6.498e-05   5.310 1.32e-05 ***
## dLoad_nonCent 2.140e+01  9.595e+00   2.231   0.0342 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.8904 on 27 degrees of freedom
## Multiple R-squared:  0.5749, Adjusted R-squared:  0.5434 
## F-statistic: 18.26 on 2 and 27 DF,  p-value: 9.654e-06

# this is what the raw data look like
plot(x = newDF$dLoad_nonCent, y = newDF$deltaMetR)

# this package plots the partial residuals
crPlot(mm7, variable = "dLoad_nonCent")

crPlot(mm7, variable = "deltaFrq2")

summary(mm7) # final model for paper

## 
## Call:
## lm(formula = deltaMetR ~ deltaFrq2 + dLoad_nonCent, data = newDF)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.17958 -0.59123 -0.06454  0.34379  3.00903 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   1.777e-01  7.507e-01   0.237   0.8146    
## deltaFrq2     3.451e-04  6.498e-05   5.310 1.32e-05 ***
## dLoad_nonCent 2.140e+01  9.595e+00   2.231   0.0342 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.8904 on 27 degrees of freedom
## Multiple R-squared:  0.5749, Adjusted R-squared:  0.5434 
## F-statistic: 18.26 on 2 and 27 DF,  p-value: 9.654e-06

Takeaways for model 3:

Deltafreq^2 and deltaLoad are both associated with deltaMetR.
Holding deltafreq constant (not necessarily at 0), an increase in deltaload is associated with an increase in deltametabolicRate
Holding deltaLoad constant (not holding load constant, and again, not holding deltaload necessarily at 0), an increase in deltaFreq^2 is associated with an increase in deltaMetabolic rate.
We found no evidence to suggest that a change in deltaArcL^2 is associated with deltametabolicRate. (this is not saying that arcLength doesn’t affect metabolic rate).

Refref: Do regression for Average percent loading vs. delta metabolic rate / 1% load etc. It’s the last page of the document Susie sent (3 different response variables).

Covariates: avg % load, order, bee size.

(DeltaMetRate / avgPercLoad) ~ avgPercLoad + order_1 + size_pc1

newDF <- within(newDF, {dmr_apl = deltaMetR / avgPercLoad})

mod1 <- lm(dmr_apl  ~ avgPercLoad + order_1 + size_pc1, data = newDF)
summary(mod1)

## 
## Call:
## lm(formula = dmr_apl ~ avgPercLoad + order_1 + size_pc1, data = newDF)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.029235 -0.014142 -0.003264  0.008463  0.053034 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)          0.1404689  0.0176079   7.978 1.86e-08 ***
## avgPercLoad         -0.0015015  0.0003086  -4.865 4.80e-05 ***
## order_1loadedSecond -0.0212737  0.0073765  -2.884  0.00779 ** 
## size_pc1             0.0012171  0.0019858   0.613  0.54527    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.01965 on 26 degrees of freedom
## Multiple R-squared:  0.6107, Adjusted R-squared:  0.5657 
## F-statistic: 13.59 on 3 and 26 DF,  p-value: 1.58e-05

mod2 <- update(mod1, .~. - size_pc1)
anova(mod1, mod2) # remove size_pc1

## Analysis of Variance Table
## 
## Model 1: dmr_apl ~ avgPercLoad + order_1 + size_pc1
## Model 2: dmr_apl ~ avgPercLoad + order_1
##   Res.Df      RSS Df   Sum of Sq      F Pr(>F)
## 1     26 0.010044                             
## 2     27 0.010189 -1 -0.00014511 0.3756 0.5453

summary(mod2) # final model for paper

## 
## Call:
## lm(formula = dmr_apl ~ avgPercLoad + order_1, data = newDF)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.030989 -0.013911 -0.001532  0.009837  0.050990 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)          0.1431985  0.0168372   8.505 4.06e-09 ***
## avgPercLoad         -0.0015573  0.0002914  -5.344 1.21e-05 ***
## order_1loadedSecond -0.0204551  0.0071702  -2.853  0.00822 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.01943 on 27 degrees of freedom
## Multiple R-squared:  0.605,  Adjusted R-squared:  0.5758 
## F-statistic: 20.68 on 2 and 27 DF,  p-value: 3.578e-06

# get p-values for avgPercLoad and order
anova(mod2,  update(mod2, .~. - order_1)) # p-value for order

## Analysis of Variance Table
## 
## Model 1: dmr_apl ~ avgPercLoad + order_1
## Model 2: dmr_apl ~ avgPercLoad
##   Res.Df      RSS Df  Sum of Sq      F   Pr(>F)   
## 1     27 0.010189                                 
## 2     28 0.013260 -1 -0.0030713 8.1384 0.008216 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

anova(mod2, update(mod2, .~. - avgPercLoad)) # p-value for avg perc load

## Analysis of Variance Table
## 
## Model 1: dmr_apl ~ avgPercLoad + order_1
## Model 2: dmr_apl ~ order_1
##   Res.Df      RSS Df Sum of Sq      F    Pr(>F)    
## 1     27 0.010189                                  
## 2     28 0.020965 -1 -0.010776 28.554 1.211e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# plot(I(deltaMetR)  ~ avgPercLoad, data = newDF)
# plot(I(deltaMetR / avgPercLoad)  ~ deltaMetR, data = newDF)

car::scatterplotMatrix(newDF[, c("avgPercLoad", "deltaMetR", "size_pc1")])

car::vif(mod1)  # looks fine

## avgPercLoad     order_1    size_pc1 
##    1.114454    1.051699    1.116654

car::vif(mod2)

## avgPercLoad     order_1 
##    1.017216    1.017216

par(mfrow = c(2,2))
plot(mod2) # possible nonlinear trend in residuals

plot(mod1, which = 4)

# update model to add a non-linear term
newDF <- within(newDF, {avgPercLoad_cent = as.numeric(scale(newDF$avgPercLoad, center = TRUE, scale = FALSE))})
newDF$apl_cent2 <- with(newDF, avgPercLoad_cent^2)

m11 <- lm(dmr_apl  ~ avgPercLoad_cent  + apl_cent2 + order_1 + size_pc1, data = newDF)
summary(m11)

## 
## Call:
## lm(formula = dmr_apl ~ avgPercLoad_cent + apl_cent2 + order_1 + 
##     size_pc1, data = newDF)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.030273 -0.010355 -0.003798  0.006817  0.055061 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)          5.247e-02  5.796e-03   9.053 2.29e-09 ***
## avgPercLoad_cent    -1.553e-03  3.088e-04  -5.030 3.45e-05 ***
## apl_cent2            2.772e-05  2.280e-05   1.216  0.23544    
## order_1loadedSecond -2.364e-02  7.564e-03  -3.125  0.00446 ** 
## size_pc1             1.321e-03  1.970e-03   0.671  0.50864    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.01948 on 25 degrees of freedom
## Multiple R-squared:  0.6324, Adjusted R-squared:  0.5736 
## F-statistic: 10.75 on 4 and 25 DF,  p-value: 3.289e-05

car::vif(m11) # vif is much better with centered variables

## avgPercLoad_cent        apl_cent2          order_1         size_pc1 
##         1.135857         1.108436         1.126036         1.118753

m11a <- update(m11, .~. - size_pc1)
anova(m11, m11a) # remove size_pc1

## Analysis of Variance Table
## 
## Model 1: dmr_apl ~ avgPercLoad_cent + apl_cent2 + order_1 + size_pc1
## Model 2: dmr_apl ~ avgPercLoad_cent + apl_cent2 + order_1
##   Res.Df       RSS Df   Sum of Sq      F Pr(>F)
## 1     25 0.0094835                             
## 2     26 0.0096541 -1 -0.00017058 0.4497 0.5086

m11b <- update(m11a, .~. - apl_cent2)
anova(m11a, m11b) # drop squared term

## Analysis of Variance Table
## 
## Model 1: dmr_apl ~ avgPercLoad_cent + apl_cent2 + order_1
## Model 2: dmr_apl ~ avgPercLoad_cent + order_1
##   Res.Df       RSS Df   Sum of Sq      F Pr(>F)
## 1     26 0.0096541                             
## 2     27 0.0101893 -1 -0.00053522 1.4414 0.2407

summary(m11b)

## 
## Call:
## lm(formula = dmr_apl ~ avgPercLoad_cent + order_1, data = newDF)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.030989 -0.013911 -0.001532  0.009837  0.050990 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)          0.0549528  0.0052157  10.536 4.58e-11 ***
## avgPercLoad_cent    -0.0015573  0.0002914  -5.344 1.21e-05 ***
## order_1loadedSecond -0.0204551  0.0071702  -2.853  0.00822 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.01943 on 27 degrees of freedom
## Multiple R-squared:  0.605,  Adjusted R-squared:  0.5758 
## F-statistic: 20.68 on 2 and 27 DF,  p-value: 3.578e-06

par(mfrow = c(2,2))

plot(m11b) # residuals look much better

plot(m11b, which = 4)

summary(m11b)

## 
## Call:
## lm(formula = dmr_apl ~ avgPercLoad_cent + order_1, data = newDF)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.030989 -0.013911 -0.001532  0.009837  0.050990 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)          0.0549528  0.0052157  10.536 4.58e-11 ***
## avgPercLoad_cent    -0.0015573  0.0002914  -5.344 1.21e-05 ***
## order_1loadedSecond -0.0204551  0.0071702  -2.853  0.00822 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.01943 on 27 degrees of freedom
## Multiple R-squared:  0.605,  Adjusted R-squared:  0.5758 
## F-statistic: 20.68 on 2 and 27 DF,  p-value: 3.578e-06

## visualize model for deltametRate/avgPercLoading

par(mfrow = c(1,1))

plot(dmr_apl  ~ avgPercLoad, col = factor(order_1), data = newDF, pch = 20)

ndf <- data.frame(avgPercLoad_cent = seq(min(newDF$avgPercLoad_cent), max(newDF$avgPercLoad_cent), length.out = 200), 
                  order_1 = as.factor(rep(c("loadedFirst", "loadedSecond"), 100)))

ndf$apl_cent2 <- ndf$avgPercLoad_cent^2

ndf$avgPercLoad <- ndf$avgPercLoad_cent + mean(newDF$avgPercLoad)
predDF <- data.frame(preds = predict(mod2, newdata = ndf), ndf)

ggplot(newDF, aes(x = avgPercLoad, y = dmr_apl)) + 
     geom_point(aes(color = order_1), , shape = 17) + 
     geom_line(data = predDF, aes(x = avgPercLoad, y = preds, color = order_1)) + 
     labs(x = "Average Load (% bodymass)", 
          y = "Change in Metabolic Rate (mL CO2 / hr) / \n Average Load (% bodymass)") + 
     scale_color_viridis( name = "Order", discrete = TRUE, end = 0.8) + 
     theme(legend.position = c(0.8,0.8))

ggsave("dmr_apl.pdf", width = 5, height = 4)

(DeltaFrq^2/avgPercload) ~ avgPercLoad + order_1 + size_pc1

newDF <- within(newDF, {df2_apl = deltaFrq2 / avgPercLoad})
par(mfrow = c(1,1))
plot(df2_apl ~ avgPercLoad, newDF)

mod1_f <- lm(df2_apl  ~ avgPercLoad + order_1 + size_pc1, data = newDF)
summary(mod1_f)

## 
## Call:
## lm(formula = df2_apl ~ avgPercLoad + order_1 + size_pc1, data = newDF)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -56.153 -25.697   3.939  22.568  51.988 
## 
## Coefficients:
##                     Estimate Std. Error t value Pr(>|t|)    
## (Intercept)         221.4973    29.3034   7.559 5.04e-08 ***
## avgPercLoad          -3.0335     0.5136  -5.906 3.14e-06 ***
## order_1loadedSecond -37.1107    12.2762  -3.023  0.00557 ** 
## size_pc1             -2.6032     3.3048  -0.788  0.43800    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 32.71 on 26 degrees of freedom
## Multiple R-squared:  0.667,  Adjusted R-squared:  0.6286 
## F-statistic: 17.36 on 3 and 26 DF,  p-value: 2.15e-06

mod2_f <- update(mod1_f, .~. - size_pc1)
anova(mod1_f, mod2_f) # remove size_pc1

## Analysis of Variance Table
## 
## Model 1: df2_apl ~ avgPercLoad + order_1 + size_pc1
## Model 2: df2_apl ~ avgPercLoad + order_1
##   Res.Df   RSS Df Sum of Sq      F Pr(>F)
## 1     26 27819                           
## 2     27 28482 -1   -663.86 0.6205  0.438

summary(mod2_f) # not final model for paper

## 
## Call:
## lm(formula = df2_apl ~ avgPercLoad + order_1, data = newDF)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -54.479 -27.491   1.816  21.916  49.750 
## 
## Coefficients:
##                     Estimate Std. Error t value Pr(>|t|)    
## (Intercept)         215.6591    28.1506   7.661 3.07e-08 ***
## avgPercLoad          -2.9140     0.4873  -5.980 2.23e-06 ***
## order_1loadedSecond -38.8616    11.9881  -3.242  0.00315 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 32.48 on 27 degrees of freedom
## Multiple R-squared:  0.6591, Adjusted R-squared:  0.6338 
## F-statistic:  26.1 on 2 and 27 DF,  p-value: 4.904e-07

# get p-values for avgPercLoad 
anova(mod2_f,  update(mod2_f, .~. - order_1)) # p-value for order

## Analysis of Variance Table
## 
## Model 1: df2_apl ~ avgPercLoad + order_1
## Model 2: df2_apl ~ avgPercLoad
##   Res.Df   RSS Df Sum of Sq      F   Pr(>F)   
## 1     27 28482                                
## 2     28 39568 -1    -11086 10.508 0.003152 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

par(mfrow = c(2,2))
plot(mod2_f) # non-linearity in residuals

plot(mod2_f, which = 4)


# update model to add a non-linear term
m22 <- lm(df2_apl  ~ avgPercLoad_cent  + apl_cent2 + order_1 + size_pc1, data = newDF)
summary(m22)

## 
## Call:
## lm(formula = df2_apl ~ avgPercLoad_cent + apl_cent2 + order_1 + 
##     size_pc1, data = newDF)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -35.430 -16.698  -1.678  19.332  43.880 
## 
## Coefficients:
##                      Estimate Std. Error t value Pr(>|t|)    
## (Intercept)          35.70197    7.30757   4.886 5.01e-05 ***
## avgPercLoad_cent     -3.27920    0.38928  -8.424 9.05e-09 ***
## apl_cent2             0.13215    0.02874   4.597 0.000106 ***
## order_1loadedSecond -48.37471    9.53594  -5.073 3.09e-05 ***
## size_pc1             -2.10870    2.48324  -0.849 0.403845    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 24.56 on 25 degrees of freedom
## Multiple R-squared:  0.8196, Adjusted R-squared:  0.7907 
## F-statistic: 28.39 on 4 and 25 DF,  p-value: 5.684e-09

car::vif(m22) # vif is much better with centered variables

## avgPercLoad_cent        apl_cent2          order_1         size_pc1 
##         1.135857         1.108436         1.126036         1.118753

m22a <- update(m22, .~. - size_pc1)
anova(m22, m22a) # remove size_pc1

## Analysis of Variance Table
## 
## Model 1: df2_apl ~ avgPercLoad_cent + apl_cent2 + order_1 + size_pc1
## Model 2: df2_apl ~ avgPercLoad_cent + apl_cent2 + order_1
##   Res.Df   RSS Df Sum of Sq      F Pr(>F)
## 1     25 15074                           
## 2     26 15509 -1    -434.8 0.7211 0.4038

summary(m22a)

## 
## Call:
## lm(formula = df2_apl ~ avgPercLoad_cent + apl_cent2 + order_1, 
##     data = newDF)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -42.409 -16.168  -2.024  16.653  47.040 
## 
## Coefficients:
##                      Estimate Std. Error t value Pr(>|t|)    
## (Intercept)          36.34581    7.22905   5.028 3.12e-05 ***
## avgPercLoad_cent     -3.18454    0.37097  -8.584 4.60e-09 ***
## apl_cent2             0.13321    0.02856   4.664 8.17e-05 ***
## order_1loadedSecond -49.88054    9.31922  -5.352 1.33e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 24.42 on 26 degrees of freedom
## Multiple R-squared:  0.8144, Adjusted R-squared:  0.793 
## F-statistic: 38.02 on 3 and 26 DF,  p-value: 1.183e-09

par(mfrow = c(2,2))

plot(m22a) # residuals look much better

plot(m22a, which = 4)

summary(m22a) # final model for paper (though we could also report the non-centered values)

## 
## Call:
## lm(formula = df2_apl ~ avgPercLoad_cent + apl_cent2 + order_1, 
##     data = newDF)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -42.409 -16.168  -2.024  16.653  47.040 
## 
## Coefficients:
##                      Estimate Std. Error t value Pr(>|t|)    
## (Intercept)          36.34581    7.22905   5.028 3.12e-05 ***
## avgPercLoad_cent     -3.18454    0.37097  -8.584 4.60e-09 ***
## apl_cent2             0.13321    0.02856   4.664 8.17e-05 ***
## order_1loadedSecond -49.88054    9.31922  -5.352 1.33e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 24.42 on 26 degrees of freedom
## Multiple R-squared:  0.8144, Adjusted R-squared:  0.793 
## F-statistic: 38.02 on 3 and 26 DF,  p-value: 1.183e-09

## visualize model for deltametRate/avgPercLoading

ndf <- data.frame(avgPercLoad_cent = seq(min(newDF$avgPercLoad_cent), max(newDF$avgPercLoad_cent), length.out = 200), 
                  order_1 = as.factor(rep(c("loadedFirst", "loadedSecond"), 100)))

ndf$apl_cent2 <- ndf$avgPercLoad_cent^2

ndf$apl <- ndf$avgPercLoad_cent + mean(newDF$avgPercLoad)

predDF <- data.frame(preds = predict(m22a, newdata = ndf, se = TRUE), ndf)

par(mfrow = c(1,1))

plot(df2_apl  ~ avgPercLoad, col = factor(order_1), data = newDF, pch = 20)

ggplot(newDF, aes(x = avgPercLoad, y = df2_apl)) + 
     geom_point(aes(color = order_1), , shape = 17) + 
     geom_line(data = predDF, aes(x = apl, y = preds.fit, color = order_1)) + 
     # geom_line(data = predDF, aes(x = apl, y = preds.fit + 1.96*preds.se.fit, color = order_1), lty = 2) + 
     # geom_line(data = predDF, aes(x = apl, y = preds.fit - 1.96*preds.se.fit, color = order_1), lty = 2) + 
     labs(x = "Average Load (% bodymass)", 
          y = "Change in wingbeat freq^2 (hz^2) / \n Average Load (% bodymass)") + 
     scale_color_viridis( name = "Order", discrete = TRUE, end = 0.8) + 
     theme(legend.position = c(0.8,0.8))

ggsave("df2_apl.pdf", width = 5, height = 4)

(DeltaArcL^2/avgPercload) ~ avgPercLoad + order_1 + size_pc1

newDF <- within(newDF, {da2_apl = deltaArcL2 / avgPercLoad})

par(mfrow = c(1,1))
plot(da2_apl ~ avgPercLoad, ylab = c("deltaArcL^2 / avgPercLoad"), newDF)

mod1_a <- lm(da2_apl  ~ avgPercLoad + order_1 + size_pc1, data = newDF)
summary(mod1_a)

## 
## Call:
## lm(formula = da2_apl ~ avgPercLoad + order_1 + size_pc1, data = newDF)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -1.023e-06 -5.865e-07  7.780e-08  4.151e-07  1.346e-06 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)          3.321e-06  6.130e-07   5.418 1.12e-05 ***
## avgPercLoad         -2.018e-08  1.074e-08  -1.878   0.0717 .  
## order_1loadedSecond -6.484e-07  2.568e-07  -2.525   0.0180 *  
## size_pc1             1.764e-08  6.913e-08   0.255   0.8006    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6.843e-07 on 26 degrees of freedom
## Multiple R-squared:  0.3183, Adjusted R-squared:  0.2397 
## F-statistic: 4.047 on 3 and 26 DF,  p-value: 0.01737

mod2_a <- update(mod1_a, .~. - size_pc1)
anova(mod1_a, mod2_a) # remove size_pc1

## Analysis of Variance Table
## 
## Model 1: da2_apl ~ avgPercLoad + order_1 + size_pc1
## Model 2: da2_apl ~ avgPercLoad + order_1
##   Res.Df        RSS Df   Sum of Sq      F Pr(>F)
## 1     26 1.2174e-11                             
## 2     27 1.2204e-11 -1 -3.0501e-14 0.0651 0.8006

summary(mod2_a) # final model for paper

## 
## Call:
## lm(formula = da2_apl ~ avgPercLoad + order_1, data = newDF)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -1.033e-06 -5.437e-07  9.477e-08  4.213e-07  1.332e-06 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)          3.361e-06  5.827e-07   5.768 3.91e-06 ***
## avgPercLoad         -2.099e-08  1.009e-08  -2.081   0.0471 *  
## order_1loadedSecond -6.365e-07  2.482e-07  -2.565   0.0162 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6.723e-07 on 27 degrees of freedom
## Multiple R-squared:  0.3166, Adjusted R-squared:  0.266 
## F-statistic: 6.255 on 2 and 27 DF,  p-value: 0.005861

# get p-values for avgPercLoad and order
anova(mod2_a,  update(mod2_a, .~. - order_1)) # p-value for order

## Analysis of Variance Table
## 
## Model 1: da2_apl ~ avgPercLoad + order_1
## Model 2: da2_apl ~ avgPercLoad
##   Res.Df        RSS Df   Sum of Sq    F  Pr(>F)  
## 1     27 1.2204e-11                              
## 2     28 1.5178e-11 -1 -2.9742e-12 6.58 0.01619 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

anova(mod2_a, update(mod2_a, .~. - avgPercLoad)) # p-value for avg perc load

## Analysis of Variance Table
## 
## Model 1: da2_apl ~ avgPercLoad + order_1
## Model 2: da2_apl ~ order_1
##   Res.Df        RSS Df   Sum of Sq     F  Pr(>F)  
## 1     27 1.2204e-11                               
## 2     28 1.4161e-11 -1 -1.9567e-12 4.329 0.04708 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

par(mfrow = c(2,2))
plot(mod2_a)

plot(mod2_a, which = 4)

# visualize model

ndf <- data.frame(avgPercLoad_cent = seq(min(newDF$avgPercLoad_cent), max(newDF$avgPercLoad_cent), length.out = 200), 
                  order_1 = as.factor(rep(c("loadedFirst", "loadedSecond"), 100)))

ndf$apl_cent2 <- ndf$avgPercLoad_cent^2

ndf$avgPercLoad <- ndf$avgPercLoad_cent + mean(newDF$avgPercLoad)

predDF <- data.frame(preds = predict(mod2_a, newdata = ndf, se = TRUE), ndf)

par(mfrow = c(1,1))

plot(da2_apl  ~ avgPercLoad, col = factor(order_1), data = newDF, pch = 20)

ggplot(newDF, aes(x = avgPercLoad, y = da2_apl)) + 
     geom_point(aes(color = order_1), shape = 17) + 
     geom_line(data = predDF, aes(x = avgPercLoad, y = preds.fit, color = order_1)) + 
     # geom_line(data = predDF, aes(x = apl, y = preds.fit + 1.96*preds.se.fit, color = order_1), lty = 2) + 
     # geom_line(data = predDF, aes(x = apl, y = preds.fit - 1.96*preds.se.fit, color = order_1), lty = 2) + 
     labs(x = "Average Load (% bodymass)", 
          y = "Change in arc length^2 (radians^2) / \n Average Load (% bodymass)") + 
     scale_color_viridis( name = "Order", discrete = TRUE, end = 0.8) + 
     theme(legend.position = c(0.8,0.8))

ggsave("da2_apl.pdf", width = 5, height = 4)

Session Info

sessionInfo()

## R version 3.3.2 (2016-10-31)
## Platform: x86_64-apple-darwin13.4.0 (64-bit)
## Running under: OS X El Capitan 10.11.6
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] plyr_1.8.4         viridis_0.3.4      visreg_2.3-0      
##  [4] effects_3.1-2      sjPlot_2.2.1       influence.ME_0.9-8
##  [7] reshape2_1.4.2     MASS_7.3-45        gsheet_0.4.2      
## [10] lme4_1.1-12        Matrix_1.2-8       car_2.1-4         
## [13] ggplot2_2.2.1     
## 
## loaded via a namespace (and not attached):
##  [1] tidyr_0.6.1        splines_3.3.2      merTools_0.3.0    
##  [4] modelr_0.1.0       shiny_1.0.0        assertthat_0.1    
##  [7] stats4_3.3.2       coin_1.1-3         yaml_2.1.14       
## [10] backports_1.0.5    lattice_0.20-34    quantreg_5.29     
## [13] digest_0.6.12      minqa_1.2.4        colorspace_1.3-2  
## [16] sandwich_2.3-4     htmltools_0.3.5    httpuv_1.3.3      
## [19] psych_1.6.12       broom_0.4.1        SparseM_1.74      
## [22] haven_1.0.0        purrr_0.2.2        xtable_1.8-2      
## [25] mvtnorm_1.0-5      scales_0.4.1       stringdist_0.9.4.4
## [28] MatrixModels_0.4-1 arm_1.9-3          tibble_1.2        
## [31] mgcv_1.8-16        DT_0.2             TH.data_1.0-8     
## [34] nnet_7.3-12        lazyeval_0.2.0     pbkrtest_0.4-6    
## [37] mnormt_1.5-5       survival_2.40-1    magrittr_1.5      
## [40] mime_0.5           evaluate_0.10      nlme_3.1-130      
## [43] foreign_0.8-67     tools_3.3.2        multcomp_1.4-6    
## [46] stringr_1.1.0      munsell_0.4.3      blme_1.0-4        
## [49] grid_3.3.2         nloptr_1.0.4       htmlwidgets_0.8   
## [52] labeling_0.3       rmarkdown_1.3      gtable_0.2.0      
## [55] codetools_0.2-15   curl_2.3           abind_1.4-5       
## [58] sjstats_0.7.1      DBI_0.5-1          sjmisc_2.2.1      
## [61] R6_2.2.0           gridExtra_2.2.1    zoo_1.7-14        
## [64] knitr_1.15.1       dplyr_0.5.0        rprojroot_1.2     
## [67] readr_1.0.0        modeltools_0.2-21  stringi_1.1.2     
## [70] parallel_3.3.2     Rcpp_0.12.9        coda_0.19-1       
## [73] lmtest_0.9-34

Sys.time()

## [1] "2017-03-16 10:12:29 EDT"

BeeResp_15Feb

Callin Switzer

February 15, 2017; Update 1 March 2017; Update March 16

Summary:

Measured variables:

Calculated variables:

Install packages and read data

Calculate variables

Create a new dataframe that calculates the changes for each bee

Calculate \(\Delta\) variables

Use PCA to combine bee sizes into a single predictor

1. First we see how \(\Delta\)load affects \(\Delta\)ArcLength accounting for bee size and treatment order

model diagnostics

Model visualization

Takeaways from model 1:

See how \(\Delta\)load affects \(\Delta\)frequency, accounting for bee size and treatment order

model diagnostics

model visualization

Takeaways from model 2:

2. We check to see how \(\Delta\)frequency and \(\Delta\)ArcLength are associated with \(\Delta\)MetabolicRate, while accounting for \(\Delta\)load

model diagnostics

model visualization

Takeaways for model 3:

(DeltaMetRate / avgPercLoad) ~ avgPercLoad + order_1 + size_pc1

(DeltaFrq^2/avgPercload) ~ avgPercLoad + order_1 + size_pc1

(DeltaArcL^2/avgPercload) ~ avgPercLoad + order_1 + size_pc1

Session Info